
d=load("data1.txt")
latitude=d(:,1)*1e-5+30.74
longitude=d(:,2)*1e-5+103.92

radiusEquatorial = 6378.1370*1000
radiusPolar = 6356.7523*1000

function rad=toRadians(degrees)
rad=degrees*pi/180
endfunction

function root=radius(x,y)
root=sqrt(power(x, 2) + power(y, 2))
endfunction

function r=earthRadius(latitude)
theta = toRadians(latitude)
a=6378.1370*1000 
b=6356.7523*1000
m = radius(power(a, 2) * cos(theta), power(b, 2) * sin(theta))
n= radius(a * cos(theta), b * sin(theta))
r=m/n
endfunction

function dis=distance(lat1,long1,lat2,long2)
dis=2* earthRadius(lat1) * asin(sqrt(power(sin((toRadians(lat2-lat1))/2),2)+cos(toRadians(lat1))*cos(toRadians(lat2))*power(sin((toRadians(long2-long1))/2),2)))
endfunction

//與初始節點相差
args=[latitude longitude zeros(length(latitude),1)+latitude(1) zeros(length(longitude),1)+longitude(1)]

result=arrayfun(@distance,args(:,1),args(:,2),args(:,3),args(:,4))

//相鄰兩點相差
lati=latitude(2:37)-latitude(1:36)
longti=longitude(2:37)-longitude(1:36)

args=[lati longti zeros(36,1) zeros(36,1)]
result=arrayfun(@distance,args(:,1),args(:,2),args(:,3),args(:,4))
plot(result)

